Persistent currents in dipolar Bose-Einstein condensates confined in annular potentials 
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We consider a dipolar Bose-Einstein condensate confined in an annular potential, with all the 
dipoles being aligned along some arbitrary direction. In addition to the dipole-dipole interaction, 
we also assume a zero-range hard-core potential. We investigate the stability of the system against 
collapse, as well as the stability of persistent currents as a function of the orientation of the dipoles 
and of the strength of the hard-core interaction. 
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I. INTRODUCTION 

Ultracold gases of bosonic and fermionic atoms have 
shown numerous fascinating effects. In some relatively 
recent experiments, it has become possible to create 
Bose-Einstein condensates of dipolar atoms [L-Q, and 
also to confine ultracold polar molecules [!, Q. Dipo- 
lar gases may give rise to novel physics as compared 
to non-dipolar ones: unlike the usual hard-core poten- 
tial, which describes the effective atom-atom interac- 
tion of non-dipolar atoms, the dipole-dipole interaction is 
anisotropic, nonlocal and it may be partly attractive and 
partly repulsive. For example, even when all the dipoles 
are aligned (under, e.g., the action of an electric or mag- 
netic field, as we consider in the present study) , the sign 
and the strength of the interaction may vary spatially, 
and can be tuned externally by changing the orientation 
of the applied field 

In some other recent experiments, toroidal trap- 
ping potentials have been realized experimentally [9l— Il3jj , 
where Bose-Einstein condensates of non-dipolar atoms 
were observed to support persistent currents 
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realization of a dipolar Bose-Einstein condensate trapped 
in toroidal/annular potentials should therefore be feasible 
experimentally, as well. A condensate in such a topolog- 
ically non-trivial confining potential is an ideal system 
for the investigation of various interesting phenomena, 
including superfluid properties, nonlinear effects, etc. 

Theoretical studies of dipolar atoms trapped in 
toroidal potentials have been performed recently [14l - [l6f . 
In Ref. [lj], Abad et al. considered a three-dimensional 
trapping geometry and used the Gross-Pitaevskii mean- 
field approximation to investigate the ground-state and 
the rotational properties of a dipolar condensate for a 
fixed orientation of the dipoles, which was chosen to be 
on the plane of motion of the atoms (say, the x-y plane) . 
It was shown that if the dipoles are pointing along, e.g., 
the x axis, there is a high concentration of them around 
the two points with coordinates (x,y,z) = (0, ±R, 0), 
with R being the mean radius of the torus, since around 
these two points the dipoles are predominantly oriented 
head-to-tail and therefore the interaction is mostly at- 
tractive. As we see below, the dipolar energy has two de- 
generate minima at these two points. As a result, as the 



strength of the dipolar interaction increases - with all the 
other parameters kept fixed - the density shows first two 
peaks around them. Eventually this two-fold symmetry 
is broken and the cloud concentrates around only one of 
them. If the dipolar interaction becomes strong enough, 
the gas collapses. In the same study it was also found 
that the gas supports metastable currents. In a more re- 
cent paper [l5[ the same authors demonstrated that this 
system exhibits Josephson oscillations and macroscopic 
quantum self-trapping. Finally, in Ref. [l6j , Zollner et al. 
considered quasi-one-dimensional motion of aligned dipo- 
lar atoms along a toroidal trap, assuming a very tight 
confinement in the transverse direction, and studied the 
few-body regime with the method of numerical diagonal- 
ization of the many-body Hamiltonian. 

In the present study we consider a dipolar 
Bose-Einstein condensate confined in a (quasi-two- 
dimensional) annular trap, and assume that all the 
dipoles are aligned along some direction, forming an 
angle with their plane of motion. In addition to the 
dipolar interaction, we consider a zero-range contact 
potential with strength g. Using the mean- field, Gross- 
Pitaevskii approximation we determine the absolute 
and/or local minima of the energy of the system for 
several values of g and G. We identify a collapsed phase, 
in which the system cannot support itself, and a stable 
one. Within the stable phase we also find a subphase in 
which the gas supports metastable, persistent currents. 

Our results show that the boundaries between the dif- 
ferent phases have a sinusoidal dependence on the angle 
0, in agreement with a simple variational model that we 
present. The atomic density distribution that we find re- 
sembles that of the study by Abad et al. When the dipole 
moment is perpendicular to the plane of motion of the 
atoms the density is axially symmetric. As the angle 
is tilted, the density develops two maxima, and eventu- 
ally the cloud concentrates around only one of these two 
maxima. We also find that, as observed in Ref. [lij for 
= 0, even when the particle density is not circularly 
symmetric, the system may support persistent currents, 
provided that the atoms are not concentrated around a 
single density maximum. This behaviour is in a sense 
reminiscent of the case of purely contact interactions, 
where, if the ground-state density localizes, the system 
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FIG. 1: (Color online) Schematic illustration of the quasi-two- 
dimensional potential. A very tight confinement is assumed 
along the z axis, while on the x-y plane the atoms are assumed 
to move in an annulus of mean radius R and width ~ o,q. 
The dipoles are oriented on the x-z plane forming an angle O 
with the x axis due to the action of an external (magnetic or 
electric) field E. 



undergoes solid-body rotation and there is no metasta- 
bility. 

In what follows we first present our model in Sec. II. In 
Sec. Ill we investigate the stability of the system against 
collapse, and in Sec. IV the stability of persistent cur- 
rents. Finally, in Sec. V we present our conclusions. 



II. MODEL 

We consider a dipolar Bose-Einstein condensate at zero 
temperature, confined in an axially-symmetric potential 



V(r) = V r (r x ) + V z (z) = X -Mu?{r x 



2„2 



(1) 



where the z direction is chosen to be the symmetry axis 
of the potential, and rj_ is the position vector on the x-y 
plane. Also, lo z and u are the frequencies of the con- 
finement along the z axis and along the direction per- 
pendicular to it. We choose lo z /uj equal to 100, in or- 
der for Tiui z to be much larger than all the other energy 
scales in the problem, in which case the motion is quasi- 
two-dimensional. Finally, we consider R/clq — 4, where 
ao = [Ti/ {Mlu)] 1 / 2 is the oscillator length corresponding 
to ui, and M is the atomic mass. Under these conditions 
V(r) describes an annular potential with mean radius R 
and width w do, as shown schematically in Fig. 1. 

Turning to the interactions, we consider both the usual 
contact potential and the dipole-dipole interaction, 
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(r, r') = U 6(r r') + ^fe , (2) 



where Uq is the matrix element for zero-energy elastic 
atom-atom collisions. When the atoms have an electric 
dipole moment d, D 2 = d 2 /(47reo), where eo is the per- 
mittivity of the vacuum, while when the atoms have a 



magnetic moment fx, D 2 — /xo/z 2 /(47r), where [io is the 
permeability of the vacuum. Finally, 9 r d is the angle be- 
tween the dipole moment and the relative position vector 
r — r' of two dipoles. 

The assumption of strong confinement along the z axis 
implies that the motion of the atoms is frozen along this 
direction. This allows us to make an ansatz for the order 
parameter of the form \&(r) = Z(z)$(r±), with Z(z) = 
e -z /2a z J (71-0,2^1/4 bgjng the ground state of the harmonic 
oscillator with frequency uj z , and a z = [K/ '(Mw^] 1 ' 2 . 
With this ansatz, one may integrate the dipole-dipole 
interaction along the z direction and derive an effectively 
two-dimensional dipolar potential. After per forming the 
integration, this potential takes the form 
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,(3) 



where <fi is the angle in cylindrical polar coordinates, r_i_ = 
r±/(2a z ), and Kq and K\ are the zeroth-order and first- 
order modified Bessel functions of the second kind. 

Integrating also the contact interaction over the z di- 
rection we find that the total effective interaction is given 
by 



Mnt,2D(r±) = g<5(r_L) + Vrfj(r_i_), 



(4) 



where g = Uq f \4>o(z)\ 4 dz = Uo/(V2na z ). Therefore, 
$(rj_) satisfies the Gross-Pitaevskii-like equation 
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+ y r (r ± ) + Vd ip ( r± )+ 5 |$( r± )| i 



*(rx) - 
= M*(r±),(5) 



where fj, is the chemical potential, and 

W-l) = / V*(t± - r'x) l*(rl)| 2 dr> ± (6) 

is the effective dipolar interaction potential. In order 
to solve Eq. (0, we have employed a fourth-order split- 
step Fourier method within an imaginary-time propaga- 
tion approach 18]. The method requires the selection 
of an initial state, which is then propagated in imag- 
inary time until a (local or absolute) minimum of the 
energy is reached and numerical convergence is achieved. 
In what follows we fix the value of the dipolar length 
add = MD 2 /(3h 2 ), to the value of add/ao — 4/3, and 
vary and g. 

III. STABILITY OF THE GROUND STATE 

We first investigate the ground state of the system. 
Since the dipolar interaction may be partly attractive, 
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FIG. 2: (Color online) Phase diagram with the orientation 
of the dipoles O in the horizontal axis and the strength of 
the contact potential g in the vertical one. The three phases 
discussed in the text are presented: the stable phase where 
there exist metastable states (green circles), the stable phase 
that does not support persistent currents (blue squares), and 
the collapsed phase (red triangles). The dashed and dotted 
lines represent the fitting functions of Eqs. and pip , re- 
spectively. 



and we also allow for negative values of the contact in- 
teraction, we start with the stability of the ground state 
against collapse. For this purpose, we choose the ini- 
tial guess for $(r_L) to be proportional to exp[— (r± — 
i?) 2 /2ag], i.e., an axially symmetric Gaussian, which is 
peaked around rj_ = R, with a width determined by ag. 
For a given 0, we solve Eq. §5§ for different values of g, 
identifying the different phases, and then we plot the cor- 
responding phase diagram g = g(O), which is shown in 
Fig. 2. 

The boundary separating the collapsed from the stable 
phase can be fitted rather well by the following sinusoidal 
function 



#c,stab/ffo ~ 11.0cos(26) - 3.5, 



(7) 



with go = D 2 /(4a z ), which is represented by the dashed 
curve in Fig. 2. As expected, the minimum value of g 
necessary to prevent the collapse decreases as the orien- 
tation angle of the dipoles increases, since in this case 
the dipole-dipole interaction becomes increasingly repul- 
sive. 

In order to get a better understanding of the above 
results, it is instructive to consider a variational order 
parameter of the form \&(r) = Z(z)Q±(s±)&i(si), where 
all the three functions are assumed to be Gaussian with 
widths determined by the oscillator length a z , and the 
variational parameters d± and di (as discussed in de- 
tail in the Appendix). Here s± is the coordinate that 
corresponds to the transverse direction and is measured 
from the minimum of the annular potential. Also, s; 



is the coordinate that corresponds to the longitudinal 
direction and is measured from one of the two minima 
of the dipolar interaction. As shown in the Appendix, 
we first integrate the dipole-dipole interaction entering 
Eq. {2J along the vertical and transverse directions, as- 
suming a z <C d±. Then, if we consider two dipoles lo- 
cated at the positions si and s[ along the annulus, for 
small values of \s r \ = |s; — s',\, we obtain an effectively 
one-dimensional dipolar potential, which is given by 
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3sin 2 (S/i?)cos 2 0]. (8) 



As one can see, Vd{si, s' { ) separates in the relative coor- 
dinate s r = si — s[ and the center-of-mass coordinate 
S = (s; + s'i)/2, as in the case of quasi-one-dimensional 
motion in a toroidal trapping potential [lj|; see also 
Refs. [H [H|. In Eq. ©, the term that depends on the 
center-of-mass coordinate S is cx 1 — 3 sin 2 (S/R) cos 2 0, 
and has two degenerate local minima, at the points with 
coordinates (x,y) = (0,±R), if ^ 90°, as expected 
from the arguments presented earlier. 

Integrating also in the longitudinal direction, under the 
assumption of a Gaussian profile, we find that the energy 
of the gas is 
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In evaluating the last term we have introduced a cut- 
off length in the integration over s r set by the oscillator 
length a z in the z direction. Also, the factor 3/2 in the 
same term comes from averaging the function sin 2 (S'/i?) 
[that enters Vcm(<S)] along the longitudinal direction. 
This factor is accurate when di is comparable to R, or in 
other words when the density in the longitudinal direc- 
tion is either homogeneous, or close to homogeneous, as 
compared to the spatial variation of Vcm, which is of the 
order of R. 

Four important conclusions follow from Eq. 

(i) First of all, this problem resembles that of a non- 
dipolar system in two dimensions interacting via an effec- 
tive contact potential with strength g c g = g + CD 2 (1 — 
3 cos 2 Q/2)/a z , where C is a dimensionless constant of 
order unity. 

(ii) If the motion in the transverse direction is frozen, 
the variational length d± may be set equal to the oscil- 
lator length ao, and di is the only free parameter. The 
motion is then quasi-one-dimensional and, for small val- 
ues of di, the kinetic energy, which scales as df 2 , dom- 
inates over the interaction energy, which scales as d\ j~ . 
Thus, in this case the total energy is always bounded 
and there is no collapse. On the other hand, if the mo- 
tion in the transverse direction is not frozen (i.e., as the 
annulus becomes wider), the system is not necessarily 
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stable against collapse: in this case, the interaction en- 
ergy scales as (did±)~ 1 , and the kinetic energy scales as 
d^ 2 and d'J 2 along the longitudinal and the transverse 
directions, respectively. Therefore, the energy is not nec- 
essarily bounded and the system may collapse, under the 
conditions that we investigate below. 

(iii) Provided that the system is stable, the energy has 
a local minimum that determines the width of the cloud 
in the transverse and longitudinal directions. 

(iv) Finally, from Eq. one may also deduce the func- 
tional form of the minimum value of <?(0) that is neces- 
sary for the stability of the gas against collapse. Since 
a z is chosen to be small, the last two terms in Eq. ([9]) are 
the dominant ones. As a result, the critical value of g is 
given by the approximate expression 



9c, stab 16 , fnr\\ -i l 

~ — [3cos(20) - 1 , 

50 3 



(10) 



which is in rather good agreement with Eq. ([7]). 

IV. PERSISTENT CURRENTS 

We turn now to the second main question of the present 
study, which is the stability of persistent currents. To an- 
swer this question we choose the initial state to be pro- 



portional to 



e !m *exp[ 



-(r ± - R) 2 /2c 



which has the 



same density distribution as before, but now with m 
units of circulation. For a given value of g and O (in 
the part of the phase diagram where the ground state 
is stable), the final state may or may not preserve the 
circulation, indicating the support or not of persistent 
currents for the particular choice of parameters. These 
two situations are represented in Fig. 2, with green circles 
and blue squares, respectively, for one unit of circulation, 
m = 1. The boundary between these two phases is given 
approximately by 



5c,met/ffo«12.0cos(26)-2.7 



(11) 



and is represented in Fig. 2 by the dotted curve. In anal- 
ogy to the case of non-dipolar atoms [22j], we have ob- 
served metastable states with up to four units of circula- 
tion, m = 4 with the minimum value of g necessary for 
the stability of the currents increasing with m. 

Again, one can get some insight into this result from 
the variational model presented in the previous section. 
For a Bose-Einstein condensate with purely contact in- 
teractions that is confined in a ring potential, the critical 
value of Uq for the existence of persistent currents with 
one unit of circulation satisfies the equation uqUq / E^ n = 
3/2, where E kin = h 2 /(2MR 2 ) and np is the density of 
the gas (homogeneously distributed) [H, [24| ■ Using this 
result for the effective coupling constant g e g, one gets 
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The above formula differs from Eq. (|10[) in the last term. 
However, this term is small, and therefore the variational 
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FIG. 3: (Color online) Two-dimensional density distribution 
of the atoms, corresponding to different points (Q,g/go) of 
the phase diagram of Fig. 2. The phase of the order param- 
eter for each configuration is also shown in the insets, with 
the associated colorbar represented on the left of the figure. 
Panels (a) and (d) correspond to metastable states with one 
unit of circulation, (b) and (c) to metastable states with two 
units of circulation, and (e), (f) correspond to current-free 
states. The dashed line in (f) is drawn to guide the eye along 
the minimum of the annular potential. 



model implies that g c ,met varies sinusoidally with 0, and 
also that the two phase boundaries are close to each 
other, in agreement with the numerical results. 

In Fig. 3 we show the two-dimensional density distri- 
bution and the phase of the order parameter correspond- 
ing to different values of and g/go within the blue 
and the green regions of the phase diagram of Fig. 2. 
When the dipole moment is oriented perpendicularly to 
the plane of motion of the atoms, the dipolar interaction 
is isotropic and purely repulsive, resembling the contact 
potential, and as a result the density is axially symmet- 
ric. Panels (a) and (b) show the density for = 90° and 
g/go = —15.2 and —14.7, corresponding to metastable 
states with one and two units of circulation, respectively. 
Since the density is axially symmetric, the expectation 
value of the angular momentum per particle for these 
current-carrying states is L/N — 1 and 2, respectively. 

Panels (c), (d) and (e) show the density for = 45° 
and g/go — —1.0,-2.8 and —2.9 respectively, corre- 
sponding to states with m = 2, 1, and 0. In this case, 
the dipolar interaction is anisotropic and as a result the 
cloud concentrates around two opposite ends of the an- 
nulus, as discussed earlier Notice that the states 

shown in (c) and (d) support persistent currents, even 
though the density is not axially symmetric. As a result, 
L/N does not have an integer value and it is equal to 
L/N = 1.87, 0.56 in panels (c) and (d), respectively. 

Finally, in panel (f) we show the density for = 15° 
and g/go = 6.0. In this case there is a spontaneous 
symmetry-breaking due to the dipolar interaction and 
the density localizes in one of the two degenerate min- 
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ima of the dipolar potential, in close analogy to the case 
of toroidal trapping [3]. The system then rotates as a 
solid body and persistent currents are not stable. Indeed, 
our results indicate that the presence of only one density 
maximum excludes the possibility of mctastability. 



V. SUMMARY AND OVERVIEW 

In this study we have considered a dipolar Bose- 
Einstcin condensate that is confined in an annular poten- 
tial. The atoms are assumed to interact via the dipolar 
interaction, plus an extra short-ranged potential. The 
dipoles are assumed to be aligned along a direction that 
forms an angle with their plane of motion. Our two 
basic results are the stability of the ground state against 
collapse and the stability of persistent currents as a func- 
tion of the angle O and of the strength of the hard-core 
potential g. 

The hard-core and the dipolar interactions combine, 
giving rise to an effective interaction that is spatially de- 
pendent. Furthermore, this effective interaction depends 
on the angle of orientation of the dipoles, which can be 
tuned externally by simply tilting the applied polarizing 
field. It is well known that the strength of the contact 
potential may be tuned via the so-called Feshbach res- 
onances [25j]. In our suggested physical system the tun- 
ability becomes possible via the change of the orientation 
of the external polarizing field. In addition, in the case 
of Feshbach resonances one may only change the value of 
the scattering length, while here the effective coupling is 
spatially dependent. 

In the phase diagram g versus we determined two 
boundaries, one delimitating the region where the system 
is stable against collapse, and another one delimitating 
the region where persistent currents are stable. Both 
phase boundaries have a sinusoidal dependence on the 
angle 0, in agreement with a simple variational model 
that we have presented. As the strength of the dipolar 
interaction decreases, the variation of g(Q) decreases and 
in the limit that it vanishes, the two phase boundaries 
which we have found become horizontal straight lines, 
corresponding to the limit of purely contact interactions. 

We have seen that metastability of the flow is possible 
even for values of the angle that give rise to inho- 
mogeneous density distributions, provided that the gas 
preserves the two-fold symmetry of the dipolar interac- 
tion and is not concentrated around only one of the min- 
ima of the dipolar potential. According to our results, 
provided that the density is not homogeneous, the pres- 
ence of two density maxima is a necessary, but not a 
sufficient condition for metastability. This situation is 
somewhat analogous to the case of a Bose-Einstein con- 
densate which is trapped in a ring potential and inter- 
acts via purely contact interactions: in this case, if the 
strength of the contact potential is sufficiently large in 
magnitude and negative, there is a spontaneous symme- 
try breaking of the axial symmetry of the Hamiltonian. 



The ground-state density localizes and it is energetically 
favorable for the system to carry angular momentum via 
solid-body rotation, in which case the system does not 
support persistent currents. 

An appealing property of this system is that one may 
sweep through all the three phases by simply changing 
the orientation angle of the dipoles, even when all the 
other parameters are kept fixed. An interesting poten- 
tial application of this fact might be a "mesoscopic su- 
perfluid ring", i.e., a device that is either "normal" or 
"superfluid" , depending on the orientation of the exter- 
nal polarizing field. Therefore, this problem is interesting 
not only for purely theoretical reasons, but also in con- 
nection with applications on e.g., precision measurements 
and nanotechnology. 
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Appendix: Evaluation of the dipolar energy 

The interaction between two dipoles with dipole mo- 
ment d located at the positions r and r' is given by 

;n ^9 1 — 3 cos 2 9 r d .. . 

y dip ,3D(r,r') = D 2 |r _ r/|3 (A.l) 

where 9 r d is the angle between the relative position vector 
r — r' and d. The relationship between D and d is given 
in Sec. II. 

First we integrate over the vertical and transverse di- 
rections, assuming that the order parameter has a prod- 
uct form *(r) = Z(z)®x(sx)$l(si), with 
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(A.3) 



The effective one-dimensional potential between two 
atoms which are located at the points s; and s[ along 
the annulus is 



Wd(*i,*i)= / ^di P ,3D(r,r')|Z(z)| 2 |$ ± ( S± )| 2 

x\Z(z')\ 2 \<P ± {s' ± )\ 2 dz d Sl _ dz'ds' ± .{AA) 

Since the confinement along the z direction is very tight, 
i.e., a z <C d±, we consider for convenience the limit a z — ¥ 
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0. In this case, 



the dipolar energy Ed of the system is given by 



2ir a z 



_> S(z r ), 



(A.5) 



i.e., we get a delta function in the integration over the 
variable z r — z — z' . In this case Vdi p ,iD is given by the 
approximate expression 



2tt D z 
~A^~dJ 



^£7(3/2, 0,2§ 2 )- 



+ -(1 - 3 sin 2 a)e s2 [2S 2 K (S 2 ) - (2s 2 - l)^i(§ 2 )] 

(A.6) 

where U is the confluent hypergeometric function, s = 
s r /(2d±), with s r = Sj — s'i being the relative coordinate, 
and sin 2 a = sin 2 (S/R) cos 2 9. 

Both for small and for large values of s the term in the 
last line of Eq. (|A.6[) is the dominant one. For \si — s[\ ^> 
d± Eq. (|A.6|) reproduces the correct asymptotic behavior 
of the dipolar potential, 



D'- 



(1 — 3 sin 2 a). 



(A.7) 



Furthermore, for small values of s, Eq. (8) follows. 

Assuming that the longitudinal wavefunction is also 
Gaussian, 
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E d = J V dip M*i,^)\Msi)\ 2 \Ms'i)\ 2 dsids' t . (A.9) 



Since the dominant contribution to this integral comes 
from small values of s, we use the expansion of 
Vdi P ,iD(sz,Sj), given by Eq. JHJ to find that 
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3 sin 2 a) 



(A.10) 



In the integration over the relative coordinate s r we have 
introduced a cut-off, which we have set equal to a z . This 
cut-off comes from the fact that in the derivation of the 
potential Vdi P) iD(sz, s[) we have taken the limit of a z —> 0, 
which gives a delta function in the integration over z — z'. 
As a result, in the initial integration, where we have the 
term (x — x 1 ) 2 + (y — y') 2 + (z — z 1 ) 2 in the denominator, 
we have set z = z' . However, in the actual calculation 
a z is small but finite. To take this into account, we have 
replaced the lower limit of s r in the integration over this 
variable by the "width" of the delta function, or in other 
words with the average value of \J ((z — z') 2 ), which is 
~ a z . 
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